function [Tsompi,Qsompi] = AR_Q_main(traces,... traces:irisFetch trace struct
    T,... %periods
    starttime,...
    endtime,...
    Q)
% calculate Q of specified periods at each channel with AR/sompi method
% Josh Crozier 2019

%Minimum and maximum number of poles of the estimated ARMA filters. 
%The minimum number of poles should be approximately twice the number of 
%clear spectral peaks.
P_min = 4;
P_max = 32;
%Minimum and maximum number of zeros of the estimated ARMA filters. 
%A small number of zeroes is useful to modelize noisy signals. 
%If q is null only AR filters are calculated.
q_min = 0;
q_max = 24;

for Tind = 1:length(T)
    
    
    % bandpass filter
%     Fstophigh = 1/(seconds(T(Tind)))/32;     %Stopband Frequency
%     Fpasshigh = 1/(seconds(T(Tind)))/4;     %Passband Frequency
%     Dpass  = 4;  % Passband Ripple
%     Dstop = 20;     %Stopband Attenuation
%     Fpasslow = 2/seconds(T(Tind));     %Passband Frequency
%     Fstoplow = 16/seconds(T(Tind));     %Stopband Frequency
%     filt = designfilt('bandpassfir','Samplerate',traces(1).sampleRate,...
%         'StopbandAttenuation1',Dstop,...        
%         'StopbandFrequency1',Fstophigh,...
%         'PassbandFrequency1',Fpasshigh,...
%         'PassbandFrequency2',Fpasslow,...
%         'StopbandFrequency2',Fstoplow,...
%         'StopbandAttenuation2',Dstop,...
%         'PassbandRipple',Dpass,...
%         'DesignMethod','kaiserwin',...
%         'ScalePassband',false);
%     %filtdelay = round(mean(grpdelay(filt,linspace(Fpasshigh,Fpasslow,10),traces(1).sampleRate)));
%     filtdelay = round((filtord(filt)/2));
    Tsompi = [];
    Qsompi = [];
    pcount = [];
    zcount = [];
    for tracein = 1:length(traces)
%     try
        
        local_starttime = max([starttime(Tind),traces(tracein).startTime]);
        startind = floor(seconds(local_starttime - traces(tracein).startTime)...
            *traces(tracein).sampleRate) + 1;
        local_endtime = min([endtime(Tind),traces(tracein).startTime...
            + seconds(traces(tracein).sampleCount/traces(tracein).sampleRate)]);
        endind = ceil(seconds(local_endtime - traces(tracein).startTime)...
            *traces(tracein).sampleRate) - 1;
%         sampleCount = endind - startind + 1;
    
        ana_sig = detrend(traces(tracein).data(startind:endind)');
%         ana_sig = filter(filt,detrend(traces(tracein).data(startind:endind)'));
%         ana_sig = [ana_sig(filtdelay+1:end),zeros(1,filtdelay)];
        [Ttemp, Qtemp, pcounttemp, zcounttemp] = AR_Q(1/traces(tracein).sampleRate,...
            ana_sig,...
            P_min, P_max, q_min, q_max,...
            seconds(T(Tind)),Q(Tind));
        Tsompi = [Tsompi;Ttemp];
        Qsompi = [Qsompi;Qtemp];
        pcount = [pcount;pcounttemp];
        zcount = [zcount;zcounttemp];

%     catch ME
%         disp(ME.message)
%         traces(tracein)
%         tracein
%         T
%         any(~isfinite(traces(tracein).data))
%     end
        f_Q_graphic(1/traces(1).sampleRate, 1./Tsompi, Qsompi, pcount, zcount, T(Tind),Q(Tind))
    end
%     f_Q_graphic(1/traces(1).sampleRate, 1./Tsompi, Qsompi, pcount, zcount, T(Tind),Q(Tind))
end

end

